1 Podsumowanie

TBA

2 Przygotowanie srodowiska i danych

2.1 Import bibliotek

library(xlsx)
library(DT)
library(knitr)
library(dplyr)
library(tidyr)
library(janitor)
library(imputeTS)
library(lares)
library(plotly)
library(caret)
library(qgraph)
library(ggforce)

2.2 Wczytanie i wstepna analiza danych

2.2.1 Wczytanie danych

raw_data <- read.xlsx(filename, 1)
raw_data <- as_tibble(raw_data)
dim(raw_data)
## [1] 6120   81

2.2.2 Krótka analiza surowych danych

Pierwsze 500 wierszy ze zbioru:

Podstawowe statystyki dla calego zbioru:

##  Dane bedace wartoscami brakujacymi (NA):  88 % 
##  Wartosc minimalna:  -1 
##  Wartosc maksymalna:  50000 
##  Liczba pacjentów:  375 
##  Okres zbierania danych:  2020-01-10 15:52:19  -  2020-01-23 09:09:23

Podstawowe statystyki dla poszczegolnych atrybutow:

2.3 Transformacja danych

2.3.1 Wstepne czyszczenie danych

Wstepne czyszczenie danych:

  • uzupelnienie kolumny PATIENT_ID
  • usuniecie pustych wierszy i kolumn
  • zmiana nazw kolumn
#raw_data[raw_data==-1]<-NA

#filling PATIENT_ID
id_filled <- raw_data %>% fill(PATIENT_ID)

#remove rows where all variables are empty
vars <- colnames(id_filled)[-(1:7)]
no_empty_rows<- id_filled[rowSums(is.na(id_filled[vars])) != length(vars), ]
no_empty_cols <- no_empty_rows[colSums(!is.na(no_empty_rows)) > 0]

#renaming columns
colnames_cleaned <- no_empty_cols %>% clean_names()

2.3.2 Brakujace wartosci

Eliminacja brakujących wartości na poziomie pacjenta obejmowała:

  • interpolację, jeżeli w kolumnie były co najmniej dwie wartości niebędące NA
  • stała wartość, jezeli w kolumnie byla dokladnie jedna wartość niebędąca NA

Jeżeli żadne z powyższych rozwiązań nie było możliwe, wartości NA zostawiono.

clean_NA<-function(column){
  not_NA_count<-sum(!is.na(column))
  if (not_NA_count>=2){ #interpolate
    column <- na_interpolation(column, option = "linear")
    column
  }

  else if (not_NA_count==1){ #constant value
    val <- first(na.omit(column))
    column[is.na(column)] <- val
    column
  }#default: leave NA values
  column
}

#for each patient:
# for each column:
#  clean_NA
cleaned<- colnames_cleaned%>% group_by(patient_id) %>% mutate_each(list(clean_NA))

#extract columns with attibutes only
attributes<-cleaned[-(1:7)]

3 Wyczyszczone dane - podsumowanie

Tabela pokazująca wyczyszczone dane:

DT::datatable(cleaned, style="bootstrap", filter = "top", rownames = FALSE, extensions = "Buttons", options = list(scrollX = TRUE, dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel')))

Podsumowanie zbioru:

all_dim <-dim(cleaned)
tmp <- cleaned %>% select(patient_id, outcome, gender) %>% group_by(patient_id) %>% summarise(outcome_count = mean(outcome), gender_count =  mean(gender))
## `summarise()` ungrouping output (override with `.groups` argument)
patients_count <- length(unique(raw_data$PATIENT_ID))-1
measurements_count <- all_dim[1]
mean_measures <- measurements_count/patients_count
genders <- tmp %>% count(gender_count) #Male 224, Female 151 # WYKRES
outcomes <- tmp %>% count(outcome_count) #201 recovered, 174 died # WYKRES
columns_count <- all_dim[2]
vars_count <- all_dim[2]-7
na_left <- round(100*mean(is.na(cleaned)), 0)

titles <- c("Liczba pacjentów", "Liczba pomiarów", "Średnia liczba pomiarów", "K", "M", "Smierc", "Zycie :(", "Liczba wierszy", "Liczba zmiennych", "Brakujace wartosci [%]")
values <- c(patients_count,measurements_count,mean_measures,1, 2, 0,1, columns_count,vars_count, na_left)
info_table <- data_frame(
  titles,
  values
)
## Warning: `data_frame()` is deprecated as of tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
knitr::kable(info_table)
titles values
Liczba pacjentów 375.00000
Liczba pomiarów 6106.00000
Srednia liczba pomiarów 16.28267
K 1.00000
M 2.00000
Smierc 0.00000
Zycie :( 1.00000
Liczba wierszy 81.00000
Liczba zmiennych 74.00000
Brakujace wartosci [%] 7.00000 
cleaned %>% distr(outcome, gender)

#wykres death_bothgenders, alive_bothgenders
invisible(remove(tmp))

plot_data <- cleaned%>%select(patient_id, gender, admission_time, discharge_time, outcome) %>%group_by(patient_id) %>% summarise_all(funs(first))
admission_plot <- ggplot() + 
    coord_cartesian() +
    scale_color_hue() +
    facet_wrap(~outcome) +
    layer(data=plot_data, 
          mapping=aes(
              x=admission_time, 
              y=discharge_time, colour=factor(gender)
              ), 
          stat="identity", 
          geom="point", 
          position=
              position_jitter()
    )
ggplotly(admission_plot)

3.1 Analiza wartosci atrybutow

Szczegółową analizę wartości atrybutów.

#data.frame(unclass(summary(tmp)), check.names = FALSE, stringsAsFactors = FALSE)
res <- attributes%>% sapply(function(x){c(min=min(x,na.rm = TRUE), median=median(x,na.rm = TRUE), mean = mean(x,na.rm = TRUE),max=max(x,na.rm = TRUE), variance=var(x,na.rm = TRUE), standard_deviation=sd(x,na.rm = TRUE), na_count = sum(is.na(x)))}) %>%  data.frame()
rd<- dim(res)
rows<-rd[1]
cols<-rd[2]


DT::datatable(res, style="bootstrap", filter = "top", rownames = FALSE, extensions = "Buttons", options = list(scrollX = TRUE, dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel')))



page_plot<-function(page_no, prow=3, pcol=3){
  print(
  ggplot(gather(attributes), na.action="na.omit", aes(value)) +   ggtitle(" title")+
    geom_histogram(bins = 10) + 
    facet_wrap_paginate
  (~key, ncol = pcol, nrow = prow, scales = "free_x", page = page_no)) 
  
}
plot_row = 4
plot_col = 4
pages<- ceiling(cols/(plot_row*plot_col))
for (i in 1:pages){
page_plot(i, plot_row, plot_col)
  
}

3.2 Korelacja miedzy danymi

“Sekcję sprawdzającą korelacje między zmiennymi; sekcja ta powinna zawierać jakąś formę graficznej prezentacji korelacji.”

plot_data <- corr_cross(attributes, 
max_pvalue = 0.05,
top = 20,
rm.na=TRUE)
ggplotly(plot_data)
names<-colnames(attributes)
Q <- qgraph(cor(attributes, use="complete.obs"), legend=TRUE, nodeNames=names, legend.cex=1.3, title ="Korelacja pomiedzy atrybutami", title.cex=8)
Q <- qgraph(cor(attributes, use="complete.obs"), layout = "spring", legend=TRUE, nodeNames=names, legend.cex=1.3, title ="Korelacja pomiedzy atrybutami",title.cex=8)

3.3 Zmiana [atrybutu/ow] w czasie

Interaktywny wykres lub animację prezentującą zmianę wybranych atrybutów w czasie.

timeline_plot <- ggplot() + 
    coord_cartesian() +
    scale_color_hue() +
    layer(data=cleaned, 
          mapping=aes(
              x=re_date, 
              y=hypersensitive_cardiac_troponin_i, group=patient_id
              ), 
          stat="identity", 
          geom="point", 
          position=
              position_jitter()
    )
ggplotly(timeline_plot)
timeline_plot <- ggplot(cleaned, aes(x=re_date, y=serum_chloride, colour=factor(patient_id), group=patient_id))  + geom_line() + geom_point() + facet_wrap(~outcome)
ggplotly(timeline_plot)
#timeline_data  <-cleaned%>%select(patient_id,re_date, hemoglobin, glucose) %>% transmute(id=patient_id, re_date=as.Date(re_date), hemoglobin, glucose) %>%group_by(id, re_date) %>% summarise(avg_hemoglobin = mean(hemoglobin), avg_glucose=mean(glucose)) %>%group_by(id) %>% mutate(day_no=re_date-min(re_date)+1, avg_hemoglobin=round(avg_hemoglobin,2), avg_glucose=round(avg_glucose,2))%>%ungroup() %>%select(day_no,avg_hemoglobin, avg_glucose)
timeline_data  <-cleaned%>%select(patient_id,re_date, hemoglobin, glucose, outcome, gender) %>% transmute(id=patient_id, re_date=as.Date(re_date), hemoglobin, glucose, outcome, gender) %>%group_by(id, re_date) %>% summarise(avg_hemoglobin = mean(hemoglobin), avg_glucose=mean(glucose), outcome=first(outcome), gender=first(gender)) %>%group_by(id) %>% mutate(day_no=re_date-min(re_date)+1, avg_hemoglobin=round(avg_hemoglobin,2), avg_glucose=round(avg_glucose,2))%>%ungroup() %>%select(day_no,avg_hemoglobin, avg_glucose, outcome, gender)

timeline_plot <- ggplot(timeline_data,aes(day_no,avg_hemoglobin, shape=factor(gender)))+geom_point(color="blue")+geom_point(aes(day_no,avg_glucose, shape=factor(gender)),color="black") + facet_wrap(~outcome)
ggplotly(timeline_plot)

3.4 Klasyfikator

klasyfikator przewidujący czy dany pacjent przeżyje (w tej sekcji należy wykorzystać wiedzę z pozostałych punktów oraz wykonać dodatkowe czynności, które mogą poprawić trafność predykcji); dobór parametrów modelu oraz oszacowanie jego skuteczności powinny zostać wykonane za pomocą techniki podziału zbioru na dane uczące, walidujące i testowe; trafność klasyfikacji powinna zostać oszacowana na podstawie kliku wybranych (i uzasadnionych) miar oceny klasyfikacji.

ml_data<- cleaned%>%group_by(patient_id)%>%summarise_all(funs(last))
ml_data<-na_mean(ml_data) #TO ZROBIC OSOBNO DLA ZBIORU TESTUJACEGO I UCZACEGO!
#rozwazyc: usuniecie tych kolumn (wierszy), w których jest dużo wartosci NA (np. powyżej 40%?)

#ml_data$outcome=as.factor(ml_data$outcome)

ml_data$outcome=factor(ml_data$outcome, 
                        labels = make.names(c("negative", "positive")))


inTraining <- 
    createDataPartition(
        # atrybut do stratyfikacji
        y = ml_data$outcome,
        # procent w zbiorze uczącym
        p = .75,
        # chcemy indeksy a nie listę
        list = FALSE)
training <- ml_data[ inTraining,]
testing <- ml_data[ -inTraining,]
ctrl <- trainControl(
    method = "repeatedcv",
    number = 2,
    repeats = 5)


fit <- train(outcome ~ .,
             data = training,
             method = "rf",
             trControl = ctrl,
             ntree = 10)
rfClasses <- predict(fit, newdata = testing)
confusionMatrix(data = rfClasses, testing$outcome)
ml_data$outcome=factor(ml_data$outcome, 
                        labels = make.names(levels(ml_data$outcome)))

rfGrid <- expand.grid(mtry = 10:30)
gridCtrl <- trainControl(
    method = "repeatedcv",
    summaryFunction = twoClassSummary,
    classProbs = TRUE,
    number = 2,
    repeats = 5)


fitTune <- train(outcome ~ .,
             data = training,
             method = "rf",
             metric = "ROC",
             preProc = c("center", "scale"),
             trControl = gridCtrl,
             tuneGrid = rfGrid,
             ntree = 30)
rfTuneClasses <- predict(fitTune,
                         newdata = testing)
confusionMatrix(data = rfTuneClasses, 
                testing$outcome)
ggplot(fitTune) + theme_bw()

3.4.1 Analizę ważności atrybutów najlepszego znalezionego modelu

3.5 Dodatkowa analiza

analiza typowa dla danych klinicznych, np.:

  • regresja logistyczna wraz z wzięciem pod uwagę czynników zakłócających (ang. confounding factors)
  • regresja Coxa (ang. Cox Proportional-Hazards Model).